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Abstract 

The human cochlea is a remarkable device, able to discern extremely small ampli- 
tude sound pressure waves, and discriminate between very close frequencies. Sim- 
ulation of the cochlea is computationally challenging due to its complex geometry, 
intricate construction and small physical size. We have developed, and are con- 
tinuing to refine, a detailed three-dimensional computational model based on an 
accurate cochlear geometry obtained from physical measurements. In the model, 
the immersed boundary method is used to calculate the fluid-structure interactions 
produced in response to incoming sound waves. The model includes a detailed and 
realistic description of the various elastic structures present. 

In this paper, we describe the computational model and its performance on the 
latest generation of shared memory servers from Hewlett Packard. Using compiler 
generated threads and OpenMP directives, we have achieved a high degree of par- 
allelism in the executable, which has made possible several large scale numerical 
simulation experiments that study the interesting features of the cochlear system. 
We show several results from these simulations, reproducing some of the basic known 
characteristics of cochlear mechanics. 
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1 Introduction 



The cochlea is the part of the inner ear where acoustic signals are transformed 
into neural pulses which are then signaled to the brain. It is a small organ, 
the size of a child's marble, which is a miracle of nature's engineering. It is 
sensitive to signals ranging over more than seven orders of magnitude, from 
a whisper to an explosion. We can hear sounds ranging in frequency from 
approximately 10 Hz to about 20 kHz, and distinguish between tones whose 
frequencies differ by less than a percent. Trained musicians, for example, are 
capable of differentiating between pure tones of 1000 Hz and 1001 Hz. 

These remarkable signal processing capabilities are achieved by a compli- 
cated mechanism involving both interaction of elastic material components 
immersed in fluid and electro-chemical processes. After decades of research a 
fascinating picture of the cochlear mechanics has emerged, but the precise na- 
ture of the mechanisms responsible for the extreme sensitivity, sharp frequency 
selectivity and the wide dynamic range of the cochlea remains unknown. The 
goal of our project is to build computational modeling tools that we hope 
will assist in the understanding of the cochlea works. We have constructed a 
comprehensive three-dimensional computational model of the fluid-structure 
interactions of the cochlea using the immersed boundary method. This is the 
first model that incorporates the curved cochlear anatomy based on real physi- 
cal measurements, uses the non-linear Navier-Stokes equations of viscous fluid 
dynamics and includes a detailed and realistic description of the various elastic 
structures present. For example, the helicoidal basilar membrane is modeled 
as an elastic shell using partial differential equations. 

We have developed a suite of software codes to support our studies of the 
cochlea using the immersed boundary method. These include codes for the 
generation of simulation input models (implemented in C++), the main im- 
mersed boundary numerical solver, and those for the analysis and visualization 
of results (implemented in Java and C++). The main workhorse in this suite 
is the general purpose immersed boundary code, which is written in C and 
requires extensive computing resources (CPU power, memory, and available 
disk space for storage of the results files). 

Our present work builds on the first author's implementation of the immersed 
boundary method for elastic shells [12]. In this paper we outline the con- 
struction of the cochlea model and present preliminary results of several large 
scale numerical experiments. These experiments reproduce some of the ba- 
sic features of cochlear mechanics and demonstrate the promise of large scale 
computational modeling approach to answering important questions about 
nature's miraculous engineering design of the cochlea. 
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The rest of the paper is organized as follows. Section 2 is a short introduction 
to cochlear mechanics, including a brief survey of some of the previous work on 
cochlear modeling. For a comprehensive up to date introduction to the subject 
see Geisler's excellent book [11]. Section 3 describes the general methodology of 
the immersed boundary method. We write down the fluid-structure interaction 
equations and outline the general first-order numerical method to solve these 
equations. The method has been extensively tested for different models of the 
elastic immersed boundary. For more detailed presentations see [31], [30] and 
[12]. We also discuss briefly the implementation of the immersed boundary 
algorithm and its optimization. Section 4 describes the construction of the 
cochlea model. In section 5 we present results of several numerical simulations. 
Our cochlea model is a work in progress and we conclude with an outline of 
future directions for this project. 



2 The Cochlea 

2. 1 Cochlear Mechanics 

The cochlea is a small snail-shell-like cavity in the temporal bone, which has 
two openings, the oval window and the round window. The cavity is filled 
with fluid and is sealed by two elastic membranes that cover the windows. 
The spiral canal of the cochlea is divided lengthwise into two passages by 
the cochlear partition consisting of the bony shelf and the so-called basilar 
membrane. These passages are the scala vestibuli and scala tympani, and 
they connect with each other at the apex, called the helicotrema. External 
sounds set the ear drum in motion, which is conveyed to the inner ear by the 
ossicles, three small bones of the middle ear, the malleus, incus and stapes. 
The ossicles function as an impedance matching device, focusing the energy 
of the ear drum on the oval window of the cochlea. The piston-like motion 
of the stapes against the oval window displaces the fluid of the cochlea, so 
generating traveling waves that propagate along the basilar membrane. 

Much of what we know about the waves in the cochlea was discovered in 
the 1940s by Georg von Bekesy [45], who carried out experiments in cochleae 
extracted from human cadavers. Von Bekesy studied the cochlea as a passive 
mechanical filter that utilizes a system of elastic components immersed in 
a fluid for analysis of incoming sounds. He observed that a pure tone input 
generates a traveling wave propagating along the basilar membrane. The wave 
amplitude rises gradually, reaching a peak at a specific location along the 
membrane, after which it decays rapidly. The peak location depends on the 
frequency of the input tone, with high frequencies peaking close to the stapes, 
and the lower frequencies further towards the apex. This "place principle" 
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is a crucial mechanism of frequency analysis in the cochlea. Resting on the 
basilar membrane is the microscopic organ of Corti, a complicated structure 
containing sensory receptors called hair cells. The hair cells detect fluid motion 
and provide input to the afferent nerve fibers that send action potentials to 
the brain. Thus, a pure tone input activates only a specific group of hair 
cells depending on its frequency, with the characteristic frequency locations 
monotonically decreasing along the basilar membrane from the stapes to the 
helicotrema. 

Von Bekesy found that the basilar membrane's elastic properties play an im- 
portant role in the wave mechanics of the cochlea. Despite its name, the basilar 
membrane is in fact an elastic shell whose compliance increases exponentially 
with length. (Unlike an elastic membrane, an elastic shell is not under inner 
tension, i.e. when it is cut the edges do not pull apart. The compliance of an 
elastic shell is defined as the amount of volume displaced per unit pressure 
difference across the surface of the shell.) Von Bekesy 's experiments further 
revealed that the wave propagates in the basilar membrane in the direction of 
increasing compliance regardless of the location of the source of excitation in 
the fluid. This phenomenon came to be known as "the traveling wave paradox" 
and it is very important because a significant part of our hearing depends on 
bone conduction, where the input to the cochlea is coming not through the 
stapes, but through the vibration of the bony walls. Bone conduction is easily 
demonstrated by placing a vibrating tuning fork in contact with the forehead, 
resulting in the subject hearing the tone of the fork frequency. 

Mathematically, the macro-mechanical system of the cochlea can be described 
by the Navier-Stokes equations of incompressible fluid mechanics coupled with 
equations modeling the elastic properties of the basilar membrane and the 
membranes of the oval and the round windows. The mathematical problem 
of solving this system of partial differential equations on a three-dimensional 
domain with intricate curved geometry is very difficult. Since the displace- 
ments of the basilar membrane are extremely small (on a nanometer scale), 
the system operates in a linear regime. Analysis shows that the waves in the 
cochlea resemble shallow water waves [24]. 

While the macro-mechanics of the cochlea break up the incoming sound into 
its frequency components, it was suggested as early as 1948 that a passive me- 
chanical system alone cannot explain the extreme sensitivity and frequency 
selectivity of the cochlea; some kind of amplification is necessary [13]. In- 
deed, analysis of cochlear macro-mechanics indicates that the traveling wave 
focusing is not sufficiently sharp, with some estimates suggesting that, at its 
threshold, human hearing is about a hundred times more sensitive than what 
would be expected from a passive macro-mechanical filter of the cochlea. 

In 1967 Johnstone and Boyle [16] utilized the Mossbauer effect to carry out 
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measurements of the basilar membrane vibrations in vivo, with more than 
100-fold improvement in resolution over von Bekesy's measurements. These 
measurements revealed that in the live cochlea the peak of the wave envelope 
is, in fact, well localized [36]. Since then the fine tuning of the cochlea has 
been confirmed experimentally in many studies, e.g. [37,39,18]. This led to 
attempts to understand the mechanical and electro-chemical processes at the 
level of the hair cells and to a renewed interest in the conjecture that the 
cochlea contains an energy source and acts as a mechanical amplifier. 

Further indirect support for the active amplification hypothesis was found 
with the discovery of the existence of otoacoustic emissions [20,17]. Otoacous- 
tic emissions may be recorded during or after acoustical stimulation using 
a sensitive microphone placed close to the ear drum. Emissions can also be 
detected when electric current is applied to the cochlea [27]. Spontaneous otoa- 
coustic emissions (SOAEs) occur in humans and other species. In the severe 
cases SOAEs are the cause of objective tinnitus, a common complaint of pa- 
tients with "ringing ears". Experiments show that a SOAE can be suppressed 
when a stimulus tone is presented at a nearby frequency. An isosuppression 
tuning curve is a curve obtained by measuring the otoacoustic emission while 
varying the amplitude and the frequency of the stimulus. In both mammalian 
and nonmammalian cochleae SOAEs and the process responsible for sharp 
frequency selectivity display similar characteristics. Measurements show that 
the isosuppression tuning curve closely resembles an ordinary tuning curve [8], 
which measures responsiveness to acoustical stimulus, leading to the conjec- 
ture that the cochlear amplifier is also the source of the SOAEs. 

Presently, understanding the nature of the active mechanism in the cochlea is 
the main open problem of hearing research. The live cochlea is a remarkable 
highly non-linear filter, but its function crucially depends on the underlying 
linear filter of the passive cochlear mechanics, which is still not sufficiently well 
understood. Basic questions about the role of geometry and the elastic prop- 
erties of the basilar membrane in cochlear mechanics remain open. Answering 
such questions is important not only to our understanding of the cochlea, but 
also in solving important engineering problems, such as a cochlear transducer 
design (see [34]). 



2.2 Cochlear Models 



Extensive research in cochlear modeling has been carried out over the years. 
Because of mathematical difficulty mostly simplified one or two-dimensional 
models that sought to incorporate some aspects of cochlear mechanics have 
been proposed. 
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Early one-dimensional "transmission line" model of the cochlea [48], [33], [10] 
has assumed that the fluid pressure is constant over a cross section of the 
cochlear channel. The fluid is assumed to be incompressible and inviscid, and 
the basilar membrane is modelled as a damped, forced harmonic oscillator with 
no elastic coupling along its length. Qualitatively, this model has been shown 
to capture the basic features of the basilar membrane response. Quantitatively, 
however, it yields large discrepancies with measurement results [46]. 

Two-dimensional models by Ranke [35] and Zwislocki [47] make similar as- 
sumptions on the cochlear fluid and the basilar membrane. Ranke's model 
uses a deep water approximation, while Zwislocki used the shallow water the- 
ory in his model. These models were further developed in [40], [23], [4], [3] and 
other works. The most rigorous analysis of a two-dimensional model with fluid 
viscosity was carried out by Leveque, Peskin and Lax in [24]. Their cochlea is 
represented by a plane (i.e. a strip of infinite length and infinite depth) and 
the basilar membrane, by an infinite line dividing the plane into two halves. 
The linearized equations are reduced to a functional equation by applying the 
Fourier transform in the direction parallel to the basilar membrane and then 
solving the resulting ordinary differential equations in the normal direction. 
The functional equation derived in this way is solved analytically, and the so- 
lution is evaluated both numerically and also asymptotically (by the method 
of stationary phase). This analysis reveals that the waves in the cochlea re- 
semble shallow water waves, i.e. ripples on the surface of a pond. A distinctive 
feature of this paper is the (then speculative) consideration of negative basilar 
membrane friction, i.e., of an amplification mechanism operating within the 
cochlea. 

Other two-dimensional models incorporate more sophisticated representations 
of the basilar membrane, using, for example, elastic beam and plate theory ([6], 
[38], [19], [2], [41], [15], [7], [14]). Three-dimensional models were considered 
by Steel and Taber [42] and de Boer [9], who used asymptotic methods and 
computations, obtaining a slightly improved fit of the experimental data. Their 
experience seems to indicate that geometry may play a significant role in the 
problem. In particular, the effect of the spiral coiling of the cochlea on the 
wave dynamics remains unresolved. It has been considered by several authors 
(see [44], [25], [43], [26]). 

With the development of more powerful computers it became possible to con- 
struct more detailed computational models of the cochlea. A two-dimensional 
computational model of the cochlea was constructed by Beyer [5]. In this 
model the cochlea is a flat rectangular strip divided into two equal halves 
by a line which represents the basilar membrane. The fluid is modelled by 
the full Navier-Stokes equations with viscosity terms, but elastic coupling 
along the basilar membrane is not incorporated. Beyer has used a modifica- 
tion of Peskin's immersed boundary method, originally developed for modeling 
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the fluid dynamics of the heart [29] . Several three-dimensional computational 
models have been reported, such as Kolston's model [21,22], intended to sim- 
ulate the micro-mechanics of the cochlear partition in the linear regime (i.e. 
near the threshold of hearing), and Parthasarati, Grosh and Nuttal's [28] hy- 
brid analytical-computational model using WKB approximations and finite- 
element methods. 



3 The Immersed Boundary Method 

The primary method used in our construction of the cochlea model is the 
immersed boundary method of Peskin and McQueen. It is a general numerical 
method for modeling an elastic material immersed in a viscous incompressible 
fluid. It has proved to be particularly useful for computer simulation of various 
biofluid dynamic systems. In this section we outline the general framework of 
the immersed boundary method. So far most applications of the method have 
treated the elastic (and possibly active) biological tissue as a collection of 
elastic fibers immersed in a viscous incompressible fluid. For details of this 
formulation of the method together with references to many applications see 
[31] and [30]. The immersed boundary framework is suitable for modeling 
not only elastic fibers, but also different elastic materials having complicated 
structure. The cochlea model makes an essential use of the immersed boundary 
method for shells, as developed in [12]. The main advantage of the immersed 
boundary method is its conceptual simplicity: the viscous incompressible fluid 
is described by the Navier-Stokes equations, the geometry of the model mirrors 
the real-life curved three-dimensional cochlear geometry and models for elastic 
and active material components can be naturally integrated. 

3. 1 The Equations of the Model 

The immersed boundary method is based on a Lagrangean formulation of 
the fluid-immersed material system. The fluid is described in the standard 
cartesian coordinates on R 3 , while the immersed material is described in a 
different curvilinear coordinate system. Let p and \x denote the density and 
the viscosity of the fluid, and let u(x, t) and p(x, t) denote its velocity and 
pressure, respectively. The Navier-Stokes equations of a viscous incompressible 
fluid are: 




(1) 



V-u = 0, 



(2) 
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where F denotes the density of the body force acting on the fluid. For example, 
if the immersed material is modeled as a thin shell, then F is a singular vector 
field, which is zero everywhere, except possibly on the surface representing 
the shell. The numerical method uses a discretization of the Navier-Stokes 
equations (1) and (2) on a periodic rectangular grid. 

Let X(q, t) denote the position of the immersed material in R 3 . For a shell, 
q takes values in a domain Q C R 2 , and X(q, t) is a 1-parameter family of 
surfaces indexed by t, i.e., X(q, t) is the middle surface of the shell at time t. 
Let f (q, t) denote the force density that the immersed material applies on the 
fluid. Then 



where 5 is the Dirac delta function on R 3 . This equation merely says that the 
fluid feels the force that the immersed material exerts on it, but it is important 
in the numerical method, where it is one of the equations determining fluid- 
material interaction. The other interaction equation is the no-slip condition 
for a viscous fluid: 



The system has to be completed by specifying the force f(q, t) of the im- 
mersed material. In a complicated system such as the cochlea the immersed 
material consists of many different components: membranes, bony walls, an 
elastic shell representing the basilar membrane, and various cells of the organ 
of Corti, including outer hair cells, which may actively generate forces. For each 
such component it is necessary to specify its own computation grid and an 
algorithm to compute its force f . It is in the specification of these forces that 
models for various system components integrate into the macro-mechanical 
model. 

3.2 The Numerical Method 

We describe here a first-order immersed boundary numerical scheme, which 
is the easiest to implement and has the important advantage of modularity: 
incorporating various models of immersed elastic material is straightforward. 

Let At denote the duration of a time step. It will be convenient to denote the 
time step by the superscript. For example u ra (x) = u(x, nAt). At the beginning 




(3) 



dx 

~dt 



u(X(q,t),t) 




(4) 
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of the n-th time step X™ and u n are known. Each time step proceeds as follows. 

(1) Compute the force f n that the immersed boundary applies to the fluid. 
For simple materials, such as fibers, this is a straightforward computation 
(see [31]). For a detailed description of a shell immersed boundary force 
computation see [12]. 

(2) Use (3) to compute the external force on the fluid F n . 

(3) Compute the new fluid velocity u n+1 from the Navier Stokes equations. 

(4) Use (4) to compute the new position X n+1 of the immersed material. 

We shall now describe in detail the computations in steps 2 — 4, beginning 
with the Navier-Stokes equations. 

The fluid equations are discretized on a rectangular lattice of mesh width h. 
We will make use of the following difference operators which act on functions 
defined on this lattice: 



DMl ) = *±Md« (5) 

D-^ = ^-f-^ (6) 

L>V(x) = ^( x + he i) - <K X - he i) ^ 

D° = ( J D?,D 2 , J D 3 ) (8) 

where i = 1,2,3, and ei, e 2 , e 3 form an orthonormal basis of R 3 . 

In step 3 we use the already known u n and F n to compute u n+1 and p n+1 by 
solving the following linear system of equations: 



' u n+l _ u r 



At 



+ f: u n k D±uA = - DV+ 1 D k D k u n+1 + F» (9) 

k=l / k=l 

D°-u n+1 = (10) 



Here u k D k stands for upwind differencing: 



n n ± 

k 



u n k Dl ul > 
ulD+ ul < 



Equations (9) and (10) are linear constant coefficient difference equations and, 
therefore, can be solved efficiently with the use of the Fast Fourier Transform 
algorithm. 



We now turn to the discretization of equations (3), (4). Let us assume, for 
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simplicity, that Q C R 2 is a rectangular domain over which all of the quantities 
related to the shell are defined. We will assume that this domain is discretized 
with mesh widths A^, Aq 2 and the computational lattice for Q is the set 

Q = {(i 1 Aq 1 ,i 2 Aq 2 ) | k = 1 • • i 2 = 1 • • .ra 2 }- 
In step 2 the force F n is computed using the following equation. 



F"(x) = ]T f"(q)^(x - X"(q))Aq (11) 

qGQ 

where Aq = Aq 1 Aq 2 and 5h is a smoothed approximation to the Dirac delta 
function on R 3 described below. 

Similarly, in step 4 updating the position of the immersed material X™ +1 is 
done using the equation 



X" +1 (q) =X"(q) + At^u m+1 (x)^(x-X"(q))/ i 3 , (12) 

X 

where the summation is over the lattice x = (hi, hj, hk), where i, j and k are 
integers. 

The function 8 h which is used in (11) and (12), is defined as follows: 
where 

|(3 - 2\r\ + ^1 + 4^1 -4r 2 ) |r| < 1 
<p{r) = \ I_0(2-| r |) l<|r|<2 

2 < |r| 

For an explanation of the construction of 5h see [31]. 

3. 3 Implementation and Optimization of Immersed Boundary Computations 



Our main simulation code implements the first-order immersed boundary al- 
gorithm described above. It is written in C and has been optimized to run on 
several platforms: the Silicon Graphics Origin 2000 parallel architecture, the 
Cray T90 vector-parallel machine and the HP V-class and Superdome systems. 
There is also a version for a PC running Windows, suitable for testing small 
models. The simulation code takes as input a set of files which contain the 
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description of the geometry and the material properties of the system. These 
files are created before the start of the simulation by the model generation 
software. Once the model has been read, the simulation enters a computation 
loop over the required number of time steps, generating output files containing 
information about the state of the simulated model. 

The complexity of an immersed boundary computation is determined by the 
sizes of the fluid and immersed boundary grids and by the size of the time 
step. To prevent fluid leakage the mesh width of the material grids is taken to 
equal approximately half the fluid grid mesh width. (For more on volume con- 
servation in immersed boundary calculations see [32]). The computation of the 
material forces is relatively inexpensive in time. The bulk of the computation 
is related to the fluid equations and to the fluid- immersed material interac- 
tion. The solution of the discretized fluid equations (9), (10) uses the Fast 
Fourier Transform algorithm. The other two demanding parts of each time 
step required development of a specialized algorithm that carefully partitions 
the fluid and the material grids into portions that are distributed amongst 
the available processors. This is done in such a way as to ensure that no 
two processors operate simultaneously on the same portion of the data. Nu- 
merical consistency conditions require reducing the time step when the space 
mesh width is decreased. A convergence study of the algorithm shows that 
the change in the time step is proportional to the change in mesh width. For 
example, decreasing the mesh width by a factor of 2 necessitates decreasing 
the time step by approximately a factor of 2 as well. Thus rescaling a model 
with a 128 3 -point fluid grid to a 256 3 -point model requires approximately 8 
times as much memory and 16 times as much CPU time [12]. 

For large scale simulations, such as that of the full cochlea model, extensive 
optimization of the immersed boundary code was necessary so as to reduce 
the elapsed time to a length that allowed useful experiments to be completed 
in a manageble time. The code is instrumented with calls to system timing 
routines, this information proving invaluable during porting and tuning of 
the software on different architectures. The largest and most recent immersed 
boundary numerical experiments have been carried out on the 32 processor HP 
Superdome installed at the Center for Advanced Computing Research (CACR) 
at Caltech. The HP 9000 Superdome at CACR contains 32 RISC processors 
arranged in a cell-based hierarchical crossbar architecture, with each cell con- 
sisting of 4 cpus with 8Gb of memory and an I/O sub-system. This architecture 
supports the shared memory programming model and the code efficiency was 
achieved primarily through the use of OpenMP parallelization directives. We 
used several software tools such as the HP CXperf and the KAI Guideview 
to examine the parallel efficiency of every section in the code and to identify 
data cache and TLB misses. Figure 1 shows the excellent scaling behavior of 
the simulation as a function of the number of allocated processors. Following 
this work, the "wall-clock" time per step of the simulation is approximately 
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1.3 seconds when running on all processors of the HP Superdome. 



Scalability of the Cochlea Simulation on the 32 
processor Hewlett Packard "SuperDome" 
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Fig. 1. Execution wall-clock time per time step vs. number of processors. 



4 The Cochlea Model 



The size of the human cochlea is about 1 cm x 1 cm x 1 cm, while the human 
basilar membrane is approximately 3.5 cm long, is very thin and very narrow: 
its width grows from about 150 /xm near the stapes to approximately 560 /mi 
near the helicotrema. Since the mesh width of the basilar membrane should 
approximately equal half the mesh width of the fluid grid (see section 3), a 
fluid grid of at least 256 3 points is necessary to adequately resolve the width 
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of the basilar membrane. The geometric model of the cochlear anatomy is 
based on measurements that include the position of the center line of the 
basilar membrane, its width and the cross-sectional area of the scalae. There 
are six surfaces in this model: the basilar membrane, the spiral bony shelf, 
the tubular walls of the scala vestibuli and the scala tympani and the semi- 
elliptical walls sealing the cochlear canal and containing the oval and the 
round window membranes (see Figure 2). The basilar membrane is modeled 
as an elastic shell following the prototype tested earlier (see [12]). The oval 
window and the round window membranes are are also modeled as elastic 
shells, but unlike in the case of the basilar membrane, the compliance of each 
of these shells is constant throughout the shell. The windows are geometrically 
identical: each is modeled by a rectangular grid, all of whose points outside a 
given circle are fixed. Hence each such grid represents a circular plate within 
a rectangular piece of a bone. 

The model building programs generate the approximate cochlear geometry, 
the material properties of various surfaces and the set of parameters describ- 
ing the desired simulation tests. The six surfaces of immersed material in the 
cochlea model are partitioned into 31 computational grids comprising approx- 
imately 750,000 points in total. In addition to the basilar membrane and the 
windows' grids described above there are 28 grids modeling the bony surfaces. 
The partitioning of the bony surfaces into 28 rectangular grids was neces- 
sary to minimize the total number of material points while maintaining an 
approximately uniform distance between these points. 

The passive cochlea is essentially a linear mechanical filter (see section 2) 
and the displacements occuring within the cochlea are too small to be seen 
without magnification. In our numerical experiments the displacements of the 
basilar membrane are measured on a nano-meter scale. The immersed bound- 
ary method is particularly suitable for such a simulation since it possesses a 
subgrid resolution (the material grid mesh width of our cochlea model is about 
20/im). Our numerical algorithm however uses the fully non-linear system of 
Navier-Stokes equations rather than the linearized system because the com- 
putational cost of solving the linearized equations in the immersed boundary 
method would be essentially the same as the cost of solving the fully non-linear 
system. 

The codes which are used to analyse and visualise the cochlear geometry and 
the simulated dynamics include a C++/OpenGL tool that runs on SGI and 
a Java/Java3D tool that runs on most platforms. These tools take as input 
the vertex coordinates for all computational grids in the cochlea model from 
the result files for each time step in the simulation. The tools generate a full 
3D rendering of the model geometry. Since the displacements occuring within 
the cochlea are very small they are color-coded to reveal the dynamics of 
the system. Our graphics tools also enable us to magnify the displacements to 
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Cochlea Model 




Fig. 2. A rendering of the geometric model of the cochlea with several parts of 
the outer shell removed in order to expose the cochlear partition consisting of the 
narrow basilar membrane and the bony shelf. The round window is located directly 
below the oval window and in this picture it is partially obscured by the cochlear 
partition. 

make them visible on the screen. An essential insight into the basilar membrane 
dynamics is provided by the plot of the normal displacement of its center line 
(see Figures 4, 5). Other Java tools display various important characteristics of 
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the system dynamics, such as the response of individual points on the basilar 
membrane as a function of time. All of the graphics tools have built-in facilities 
for generating animation. 

The construction of the full cochlea model has been undertaken in stages 
with the individual components tested separately prior to the model assembly. 
The basilar membrane model is described in detail in [12]. The oval window 
functions as the input window of the cochlea. We simulate the pressure of the 
stapes against it by prescribing an external force vector field on the window 
grid. This force is orthogonal to the surface of the window and for a pure tone 
input its magnitude varies sinusoidally. No force is prescribed on the round 
window. 

Since many of the numerical experiments with the full cochlea model require 
days of computing we have made a small modification of the cochlear anatomy 
in our model. In the altered model the cochlea is packed more tightly and the 
whole structure fits in a half-cube of size 1cm x 1cm x 0.5cm rather than a 
1 cm 3 cube, making it possible to use a fluid grid of 256 x 256 x 128 points. 
This configuration requires about 1.5 Gigabytes of main memory and the total 
time needed to simulate a single time step is reduced approximately in half. 
The main difference between the "half-cube model" and our original model is 
that in the former model the basilar membrane if flat. I.e. its spiraling shape 
is completely contained within a plane. All of our test simulations reported in 
the next section have been carried out with the half-cube model. 



5 Numerical Experiments 

The good efficiency of the immersed boundary solver allowed us to complete 
several large numerical experiments. In this section we present preliminary 
results from four such experiments. We are continuing this work and a more 
detailed exposition of the collected data and its analysis will be published in 
our next paper. 

Each of the experiments presented here consisted of applying a pure tone in- 
put of a given frequency at the stapes and studying the subsequent motion of 
the basilar membrane. The input at the stapes was simulated by specifying an 
external periodic force vector field on the oval window grid. A very small time 
step of 30 nano-seconds was chosen to guarantee numerical stability and good 
detail. The choice of the time step was made as a result of the convergence 
study of the system carried out in [12]. Our initial aim in the numerical exper- 
iments was to reproduce the qualitative characteristics of cochlear mechanics 
predicted by asymptotic analysis and previously reported computational mod- 
els. Consequently, we have attempted to capture the steady state response of 
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the basilar membrane to pure tones and in each of our experiments we have 
run the system for the duration equivalent to several input cycles. For exam- 
ple, for a 10 kHz input tone we have run the system for more than 30,000 
time steps, equivalent to 0.9 msec of total simulated time. On a 16 processor 
Super Dome this computation required more than 20 hours to complete. Every 
10th time step of the simulation the program generated an output file contain- 
ing the instantaneous position of the computation grids. The total amount of 
storage required is measured in tens of Gigabytes of disk space. Correspond- 
ingly, the 2 kHz experiment required five times as much computational time 
and storage. 




Fig. 3. A close-up snapshot of the traveling wave propagating along the basilar 
membrane. The magnitude of the basilar membrane displacement has been amplified 
in the normal direction. 

We have verified that all of our experiments have been carried in the system's 
linear regime, i.e. the input force at the stapes was kept very small, resulting 
in the basilar membrane displacements on the order of nano-meters. Indeed, 
increasing the force by a factor of 10 resulted in basilar membrane displace- 
ments almost exactly ten times bigger. Since in each experiment the system 
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was started from rest, we have observed initial oscillatory transient response, 
which was followed by the smoothing of the traveling wave. A close-up snap- 
shot of the traveling wave propagating along the basilar membrane is shown 
in Figure 3. The wave magnitude has been amplified in the direction of the 
normal to the basilar membrane. 




Fig. 4. A snapshot of the center line of the basilar membrane response and wave 
envelope. The top plot shows the response to the 15 kHz input sound. The wave 
snapshot is taken after 55000 time steps (1.65 msec). The wave envelope was com- 
puted over time steps 45,000 - 55,000. The bottom plot shows a 10 kHz experiment. 
The wave is shown after 30,000 time steps and the envelope was computed over 
time steps 20,000 - 30,000. The unit of vertical scale is 10 _5 cm. 

Much information about the traveling wave is revealed in the plot of the 
centerline of the basilar membrane. Four such plots are shown in Figure 4 and 
Figure 5. Figure 4 shows the results of the experiments with input frequency of 
15 kHz (top plot) and 10 kHz (bottom plot), and Figure 5 shows the results of 
the 5 kHz and the 2 kHz experiments. In each plot an instantaneous position of 
the centerline of the basilar membrane is shown, as well as the wave envelope 
computed over a range of time steps. 

Our experiments reproduced the following characteristic features of cochlear 
mechanics: In each instance, in response to a pure tone input frequency, we 
have observed a traveling wave propagating from the stapes in the direction 
of the helicotrema. The amplitude of the wave is gradually increasing until 
it reaches a peak at a characteristic location along the basilar membrane de- 
pending on the input frequency. The speed of the wave is sharply reduced 
as it propagates further along the basilar membrane. The higher the input 
frequency, the closer the peak of the wave is to the stapes. Furthermore, after 
reaching the peak the wave drops off sharply, essentially shutting down. Since 
no active mechanism has been incorporated into the present model yet, the 
observed traveling wave is not sharply focused. 



17 



Fig. 5. A snapshot of the center line of the basilar membrane response and wave 
envelope. The top plot shows the response to the 5 kHz input sound. The wave 
snapshot is taken after 70000 time steps. The wave envelope was computed over 
time steps 60,000 - 70,000. The bottom plot shows a 2 kHz experiment. The wave 
is shown after 100,000 time steps and the envelope was computed over time steps 
100,000 - 130,000. The unit of vertical scale is 10 _6 cm and 10 _7 cm in the top and 
in the bottom plots, respectively. 

It is interesting to note that while the computed traveling wave is smooth, its 
envelope is generally not smooth. The 10 kHz experiment shown in Figure 4 
is a particularly striking example, with the envelope turning 90° sharply in 
the vertical direction approximately at the 0.75 cm location along the basilar 
membrane. We believe that further testing is necessary and more data must 
be collected from numerical experiments before further conclusions about the 
shape of the traveling wave envelope can be drawn. 

The interested reader is invited to view several animations of our results by 
visiting our web site[l]. 



6 Summary and Conclusions 

We have constructed a comprehensive three-dimensional computational model 
of the passive cochlea using the immersed boundary method. Extensive opti- 
mization and parallelization made it possible to complete several large scale 
numerical simulations on a 32-processor shared memory HP Superdome com- 
puter. Together with the previous demonstration of the traveling wave paradox 
(see [12]), the pure tone experiments reported in this paper capture the most 
important properties of the cochlear macro-mechanics. 

We would like to note that our results are of somewhat preliminary nature. 
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We are continuing to test the model and plan to complete more numerical 
experiments in order to compare our simulation results with the available ex- 
perimental and modeling data. We believe our results demonstrate the promise 
of large scale computational modeling approach to the study of cochlear me- 
chanics. 

The cochlea is a very small and delicate organ and it is very difficult to study 
experimentally. Many of the important questions of cochlear mechanics are 
mathematically very complicated, but they can be studied using numerical 
simulation. Two such general questions concern the effects of the geometry 
and of the elastic properties of various cochlear components on the dynamics. 
It is easy to test cochlear models with different geometries such as a straight- 
ened out model. Testing different elastic models for the cochlear components 
is somewhat more involved since it generally requires testing each such model 
separately within the immersed boundary framework. It is nevertheless a fea- 
sible project, which is straightforward for the expert. 

Going beyound the questions of passive cochlear macro-mechanics we would 
like to incorporate an active mechanism into our model. The function of the 
live cochlea undoubtedly depends on the combination of passive cochlear me- 
chanics with an active mechanism. A comprehensive model of the passive 
cochlea is therefore a necessary first step towards modeling the active live 
cochlea. It is interesting to note that refining the mesh width of the compu- 
tational grids by a factor of 2 yields a material mesh width of approximately 
10/im. This mesh width is small enough to allow much of the structure of 
the organ of Corti to be incorporated into the model. Such a refined model is 
not yet feasible since it would require approximately 16 times more comput- 
ing power than our present model. Indeed, the large scale of computations is 
presently the biggest obstacle to progress in cochlear modeling, but this is, no 
doubt, a temporary obstacle. Continuing progress in hardware and software 
will make a construction of even the refined cochlear model possible soon. 
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